The role of pore fluids in supershear earthquake ruptures

The intensity and damage potential of earthquakes are linked to the speed at which rupture propagates along sliding crustal faults. Most earthquakes are sub-Rayleigh, with ruptures that are slower than the surface Rayleigh waves. In supershear earthquakes, ruptures are faster than the shear waves, leading to sharp pressure concentrations and larger intensities compared with the more common sub-Rayleigh ones. Despite significant theoretical and experimental advances over the past two decades, the geological and geomechanical controls on rupture speed transitions remain poorly understood. Here we propose that pore fluids play an important role in explaining earthquake rupture speed: the pore pressure may increase sharply at the compressional front during rupture propagation, promoting shear failure ahead of the rupture front and accelerating its propagation into the supershear range. We characterize the transition from sub-Rayleigh to supershear rupture in fluid-saturated rock, and show that the proposed poroelastic weakening mechanism may be a controlling factor for intersonic earthquake ruptures.


S1. Mathematical model for dynamic ruptures in poroelastic media
The mass conservation equation for the fluid is given by: written in terms of the pore pressure, p, and volumetric strain, ε v . The model parameters are the intrinsic permeability of the porous rock, k, the fluid viscosity and density, η f and ρ f , the specific storage of the saturated rock, S ϵ = ϕχ f + (α B − ϕ)χ s , where χ f and χ s are, respectively, the fluid and solid compressibilities, and the rock porosity, ϕ. The Biot coefficient is α B . Neglecting fluid accelerations and relative fluid-solid accelerations, the internal equilibrium equation for the solid is given by: where ρ b is the bulk density of the saturated rock. The Cauchy stress tensor is decomposed as the sum of an inviscid (elastic) part plus a viscoelastic one: We adopt a Saint Venant-Kirchhoff hyperelastic material for the inviscid part whose stress stensor is written in terms of the infinitesimal strain tensor as: where G is the shear modulus and λ is the first Lamé constant, λ = Eν/(1 + ν)(1 − 2ν), where E is the Young modulus and ν the Poisson ratio. For a Kelvin-Voigt material, the viscous component of the Cauchy stress tensor reads: The total stress tensor is given by: Finally, neglecting cohesion, the fault frictional strength is given by: where σ ′ n is the effective normal compression on the fault, and µ is the friction coefficient. The effective contact pressure σ ′ n is given by σ ′ n = p − T n , with T n being the contact pressure between the fault edges (compressive pressures are positive), and p is the fault pore pressure. Its value is chosen as the maximum on both sides of the fault, p = max(p − , p + ) (Jha & Juanes, 2014). The fault remains locked when the shear stress acting on the fault is lower than the frictional strength, < τ f ; otherwise, it slips. We model flow across the fault through a simple mass flux exchange between the two contact surfaces defining the fault. Denoting by p ± the pressures on either side of the fault and by q i n ± the inward mass fluxes per unit area, we approximate the latter through an effective fault transmissibility, T f , as q i n − = T f (p + − p − ), and q i n + = T f (p − − p + ). The above jump condition couples the pore pressures on both sides of the fault, allowing a transition from essentially no-flow (T f → 0) towards pressure continuity (T f → ∞) as T f increases. In this study, we consider small transmissibility, T f = 10 −14 s/m, so that pressure diffusion across the fault is negligible for the time scales studied.
The fault constitutive behaviour is governed by the Dieterich-Ruina aging law, where the friction coefficient evolves according to: and the state variable is as follows: In the above equations, V is the slip speed, and (µ 0 , a, b, D c , V 0 ) are model parameters.
For slip-weakening faults the coefficient of friction evolves from its static value, µ s , to the dynamic one, µ d , over a characteristic distance, d 0 , as where |d| is the total accumulated slip.

S2. Simulation setup
We simulate earthquakes in poroelastic media using two configurations and triggering mechanisms: 1. For Figures 1-4 of the main text: Earthquakes triggered by fluid injection near a strike-slip fault. The injection-induced model is a 2-D horizontal domain with a plane-strain poroelastic medium and an impermeable fault with rate-and-state friction, as the one described in [1] (Supplementary Figure S1). The fault is locked at the beginning of injection, with acting shear stress and frictional strength determined by anisotropic far-field tectonic stresses, σ x >σ y . Fluid injection near the center of the fault increases the pore pressure, triggering an earthquake when the acting shear stress equals frictional strength. We control the rupture propagation speed by varying the released stress during rupture (changing the confining stresses (σ x , σ y ) and their ratio, as well as rock stiffness and poroelastic coupling (through the total system compressibility or specific storage, S ϵ , and the Biot coefficient, α B ). The main field variables we use to characterize dynamic ruptures are the pore pressures, displacements, and accelerations. We also track the evolution of pore pressures on both sides of the fault, slip speed, and accumulated slip. The latter is used to estimate rupture propagation speed (Section S3). The model parameters used in the calculations leading to Figures 2, 3, 4a, and 4b of the main text are summarized in Tables S1-S3, respectively.

For Figures 5 of the main text:
Earthquakes triggered by tectonic loading. We use this idealized scenario to investigate the validity of the Andrews' seismic ratio when computed using effective normal stresses. The tectonic model is a 2D horizontal domain that encloses a 4 km length strike-slip fault parallel to the x-axis, with a plane-strain poroelastic medium and an impermeable fault with a slip-weakening friction (Supplementary Figure S2). We generate the initial conditions of confinement (normal stress acting on the fault by imposing a displacement in the y-direction along the upper and lower boundaries, negative and positive respectively. This imposed displacement leads to an initial effective normal stress (σ ′ n,0 ≈9.15 MPa). From the initial state of normal confinement, we gradually increase the applied shear stress keeping by imposing a time dependent x-displacement along the horizontal boundaries (a shear rate), and a triangular law on the vertical boundaries (Supplementary Figure  S2). The gradual increase of shear stress acting on the fault eventually leads to the nucleation of earthquakes at a small weaker patch at the center of the fault. We trigger failure by reducting assuming a smaller static friction coefficient at a small seed region around the center of the fault. The model parameters used in the calculations leading to Figure 5 of the main text are summarized in Table S4.

S3. Estimating the speed of rupture propagation
We compute the rupture propagation speed during a simulation, V R (t), using simple first-order differencing of the front positions over time (Supplementary Figure S3 below). The position of the rupture front is calculated from the profile of accumulated slip along the fault. It is defined as the coordinate along the fault of the first node whose accumulated slip is larger than a threshold value (we take 10 −4 m). We then compare rupture front positions between time steps separated by a lag that reduces differencing noise (i.e. we compute the rupture speed at step i by comparing the front position at step i and the position at step i − 50), and divide the distance travelled by the elapsed time: where V R,i is the rupture propagation speed at time level t i , and X f,i is the position of the rupture front at that same time level (coordinate along the fault). Rupture speeds calculated using finite differences are necessarily oscillatory. To reduce high-frequency oscillations in the calculation, we filter the raw rupture speeds using a moving average (Supplementary Figure S3 below). In Figure S3 we show two examples of evolution of front speed (these correspond to two data points in Figure 4 of the main text, one sub-Rayleigh and one supershear). We show the variation of the relative speeds, V R,i /C S,0 , throughout the simulation, for E = 50 GPa and α B = 0.2 (panels a and b), and E = 50 GPa and α B = 0.9 (panels c and d), which propagate, respectively, in sub-Rayleigh and in supershear regimes. In panels (b) and (d) we highlight the average value of rupture speed that we has been considered for Figure 4,V R /C S,0 . Figure S4 shows the evolution in time of rupture speeds for all the simulation cases reported in Figure 4, where we set four values of the Young modulus, E = 35, 50, 75, and 90 GPa and several values of the Biot coefficient α B , which characterize the parametric transition from sub-Rayleigh to supershear propagation. For reference, we show here again Figure 4a of the main text, which summarizes the rupture speeds as a function of rock properties E and α B (shown here as Figure S5).

S4. Impact of grid resolution on rupture speed
To assess the impact of grid resolution on the reported behavior of the rupture speeds, we performed a refinement study for a set of simulations from those shown in Figure 4 of the main manuscript, namely those with E = 75 GPa and several values of the Biot coefficient ( Figures S6 and S7). We used computational grids with grid sizes along the fault of δs = 2, 4, and 10 m. Figure S6 compares the evolution of rupture speed for four values of the Biot coefficient and the three refinement levels. While the different grid resolutions yield slightly different evolution of the rupture speed, the asymptotic rupture speeds, which are the relevant output for our study, are quite robust. The main differences among grid refinement levels arise at the parametric transition to supershear rupture (α B = 0.55, Figure S7). Here the coarse and fine level lead to subshear rupture while the intermediate grid level yields a supershear rupture. We believe this behavior is consistent with the idea that the parametric transition is truly a bifurcation with discontinuous behavior, rather than a smooth transition to supershear as the parameter (the Biot coefficient) increases. a b c d Figure S3: Calculation of rupture propagation speeds. We estimate rupture speeds using firstorder time differencing of the front locations along the fault, as determined from the evolution of the profile of accumulated slip (panels a and c). The rupture speeds calculated using finite differences are noisy, so we estimate average speeds using moving averages (panels b and d). In this figure we show the calculation of rupture speeds for two of the data points shown in Figure  4 of the main text: those with parameters E = 50 GPa and α B = 0.2 (panels a and b), and E = 50 GPa and α B = 0.9 (panels c and d). Figure S5: Summary of rupture front speeds from the previous figure (Fig. S4). This is also shown as Figure 4a of the main text.

S5. Estimating the nominal and effective seismic ratios
To illustrate the calculation of the seismic ratio, we study the changes in frictional strength and on applied shear stress during rupture propagation in earthquakes driven by tectonic loading (Supplementary Figure S4), For our slip-weakening fault, the Andrews seismic ratio, S, is simply defined as [2,3]: where µ s and µ d are the static and dynamic friction coefficients, σ ′ n is the effective normal stress, and τ a is the applied shear rate. In order to envisage the seismic ratio, we have renamed the numerator and the denominator as S1 and S2, respectively. In Supplementary Figure S4 we compare the profiles of shear stress and frictional strength along the fault, for the same overall simulation parameters, and two contrasting cases in terms of poroelastic coupling: the purely elastic case α B =0 (panel a) and the poroelastic case with α B =0.95 (panel b). Calculating the seismic ratio for the elastic case is immediate ( Supplementary  Figure 4a), and the sub-Rayleigh to supershear transition (Supplementary Figure 6) is associated to a seismic ratio of S=1.77 (as noticed by Andrews [2]). The effective seismic ratio, S ′ , is defined using effective quantities at the front in Supplementary Figure 4b